Medical imaging is one of the most useful diagnostic tools available to medicine. The invention presented here relates to one of the most popular imaging techniques belonging to the emission tomography category: positron emission tomography (PET). This medical imaging technique allows us to look inside a person and obtain images that illustrate various biological processes and functions. In this technique, a patient is initially injected with a radiotracer, which contains bio-chemical molecules. These molecules are tagged with a positron emitting radioisotope, and can participate in physiological processes in the body. After the decay of these radioisotope molecules, positrons are emitted from the various tissues of the body which have absorbed the molecules. As a consequence of the annihilation of the positrons, pairs of gamma photons are produced and are released in opposite directions. In PET scanners, these pairs of photons are registered by detectors and counted. A pair of detectors detecting a pair of gamma photons at the same time constitutes a line of response (LOR). A count of photons registered on a certain LOR will be called a projection. Data associated with annihilation events along different LORs are collected and processed. A given set of projections is mostly formed as a so-called sinogram based on their corresponding LORs.
The goal of the PET is to reconstruct the distribution of the radiotracer in the tissues of the investigated cross-sections of the body based on a set of projections from various LORs obtained by the PET scanner. The problem formulated in this way is called an image reconstruction from projections problem and is solved using various reconstruction methods. Because of the relatively small number of annihilations observed in a single LOR, the statistical nature of the measurements performed has a strong influence and must be taken into account. Recently, some new concepts regarding reconstruction algorithms have been applied to emission tomography techniques (i.e. to positron emission tomography (PET)), with statistical approaches to image reconstruction being particularly preferred (see e.g. [1], [2]). The standard reconstruction method used in PET is the maximum likelihood—expectation maximization (ML-EM) algorithm, as described for example in [3][4]. In this algorithm an iterative procedure is used in the reconstruction process, as follows:
                                          f            l                          t              +              1                                =                                    f              l              t                        ⁢                          1                                                ∑                  k                                ⁢                                  a                  kl                                                      ⁢                                          ∑                k                            ⁢                                                a                  kl                                ⁢                                                      λ                    k                                                                              ∑                                              l                        ′                                                              ⁢                                                                  a                                                  kl                          ′                                                                    ⁢                                              f                                                  l                          ′                                                t                                                                                                                                ,                            (        1        )            where: ƒl is an estimate of the image representing the distribution of the radiotracer in the body; l=1, . . . , L is an index of pixels; t is an iteration index; λk is the number of annihilation events detected along the k-th LOR; akl is an element of the system matrix.
Modification of this method, i.e. using ordered subset expectation maximization (OSEM), and improvements in computing speed, have allowed iterative algorithms to be used in the standard clinical practice of PET devices.
This algorithm is presented in the literature as being more robust and flexible than analytical inversion methods because it allows for accurate modeling of the statistics of the measurements obtained in the PET, i.e. the Poisson statistical distribution of the annihilations detected on the LORs.
The image processing methodology used in this algorithm is consistent with the algebraic image reconstruction scheme, where the reconstructed image is conceptually divided into homogeneous blocks representing pixels. In this algebraic conception, the elements of the system matrix αkl are determined for every pixel l separately, for every annihilation event λk detected along the k-th LOR. Bearing in mind the non-zero width of the radiation detector, it is easy to ascertain the set of image blocks that have an influence on the formation of the measurement λk. As can be seen from FIG. 2.a, as the ray passes through the test object (the image), all the squares through which part of the ray passes are taken into consideration. The next step is to consider the contribution made by each image block to the way in which each LOR passes through. This contribution can be calculated, for example, by counting the sub-blocks of a given block l through which the LOR k passes. Unfortunately, algebraic reconstruction problems are formulated using matrices with very large dimensionality. An algebraic reconstruction algorithms are thus much more complex than analytical methods.
The author of this invention has devised a new statistical approach to the image reconstruction problem, which is consistent with the analytical methodology of image processing during the reconstruction process. The problem as formulated by the author of this invention can be defined as an approximate discrete 2D reconstruction problem (see e.g. [5]). It takes into consideration a form of the interpolation function used in back-projection operations. The preliminary conception of this kind of image reconstruction from projections strategy for transmission tomography, i.e. x-ray computed tomography (CT), is represented in the literature only in the original works published by the author of this invention, for parallel scanner geometry (see e.g. [6]), for fan-beam geometry (see e.g. [5]) and for spiral cone-beam tomography (see e.g. [7]). The reconstruction procedure used for a spiral cone-beam CT scanner has been patented in 2016 (see U.S. Pat. No. 9,508,164 B2). Thanks to the analytical origins of the reconstruction method proposed in the above papers, most of the above-mentioned difficulties connected with using algebraic methodology can be avoided. Although the proposed reconstruction method has to establish certain coefficients, these can be pre-calculated and, because of the small memory requirements, can be stored in memory. Generally, in algebraic methods, the coefficients aid are calculated dynamically during the reconstruction process, because of the huge dimensionality of the matrix containing these elements of the system. The analytical reconstruction problem is formulated as a shift-invariant system, which allows the application of a FFT algorithm during the most demanding calculations, and in consequence, significantly accelerates the image reconstruction process.
Measurements obtained from CT are subject to statistics consistent with the Poisson distribution, but are then transformed by a log function, and so in transmission tomography, the preferred reconstruction approaches, formulated according to the ML method, are based on a weighted least squares estimate. However, in the case of emission tomography, e.g. the PET, the measurements obtained from the scanner are subject directly to statistics consistent with the Poisson distribution. This means that the preferred approaches for this imaging technique (using the ML method) are based on the Kullback-Leibrer divergence, and the EM algorithm associated with it. This invention is strictly concerned with the analytical reconstruction approach previously devised for the transmission CT technique, and the adoption of this solution to the emission PET imaging technique, using an EM algorithm with an analytical scheme of image processing. Previously, an attempt has been made to develop the idea of an iterative statistical approach with an analytical image processing framework for the CT technique in the direction of the EM method [8]. However, that idea is not consistent with the EM method (despite the authors' intentions), which implies that the reconstruction method proposed there is not optimal in so far as the form of the statistical conditions present in the PET imaging technique is concerned. Moreover, the iterative reconstruction proposed in that paper does not exploit the possibilities of FFT algorithms, which makes the method proposed in that paper computationally highly inefficient.